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We report here, for the first time, an observed instability of a Kelvin-Helmholtz (KH) nature occurring 
in a fully three-dimensional (3D) current- vortex sheet at the fan plane of a 3D magnetic null point. The 
current-vortex layer forms self-consistently in response to foot point driving around the spine lines of the 
null. The layer first becomes unstable at an intermediate distance from the null point, with the instability 
being characterized by a rippling of the fan surface and a filamentation of the current density and vorticity 
in the shear layer. Owing to the 3D geometry of the shear layer, a branching of the current filaments and 
vortices is observed. The instability results in a mixing of plasma between the two topologically distinct 
regions of magnetic flux on either side of the fan separatrix surface, as flux is reconnected across this surface. 
We make a preliminary investigation of the scaling of the system with the dissipation parameters. Our results 
indicate that the fan plane separatrix surface is an ideal candidate for the formation of current- vortex sheets 
in complex magnetic fields and, therefore, the enhanced heating and connectivity change associated with the 
instabilities of such layers. 
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I. INTRODUCTION 

In recent years three-dimensional (3D) magnetic null 
points (points in space where the magnetic field strength 
is zero) have become increasingly recognized as im- 
portant sites for energy release and magnetic topology 
change through magnetic reconnection. They are found 
in abundance in the lower solar atmosphere^ whilst dur- 
ing active times of the solar cycle they have been inferred 
to be involved in solar jets^, flux emergence^, flare 
brightening^ and magnetic breakoutP. In the Earth's 
magnetosphere their existence has been confirmed in the 
reconnecting magnetotail current sheet through in situ 
measurements 8 whilst clusters of nulls have also been 
seen in the polar cusp regions during global magneto- 
spheric simulations^. 

The magnetic field near to a 3D null can be character- 
ized by two main topological structures. The fan plane: 
a surface of field lines which emanate from (approach to- 
ward) the null and the spine lines: two field lines which 
approach toward (recede away from) the null. Reconnec- 
tion at 3D nulls can take place in several ways. Twisting 
motions about the spine/fan create a current sheet along 
the fan/spine leading to a rotational slippage known as 
torsional fan/spine reconnection, respectively!^. Shear- 
ing motions of the spine or fan in an incompressible 
plasma give rise to cur rent accumulation at the fan and 
spine respective ^ 11 1 12 1 known as fan and spine reconnec- 
tion. However, when the incompressibility assumption is 
relaxed, the plasma pressure gradient is too weak within 
the planar (spine or fan) current sheet to oppose the col- 
lapse of the field around the null (driven by the Lorentz 
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force), and a current sheet forms at an angle to the 
spine and fan within which spine-fan reconnection takes 
plac e 1 1Q l 13 -ifl5] Lastly, it has recently been suggested^ on 
the basis of an infinite-time singularity that a new regime 
exists for initially asymmetric nulls where shear builds 
aligned to the minor fan plane axis, resulting in a cur- 
rent sheet with a strong component parallel to the spine. 

In the cases of externally driven fan and torsional fan 
reconnection the current layer that forms at the fan sep- 
aratrix surface includes both a sheared magnetic and ve- 
locity field component. Studies of both regimes have un- 
til now focused purely on the formation of these smooth 
current layers. However, such a shear layer configuration 
(known as a current- vortex sheet) is know to be unstable 
to shear flow and resistive instabilities leading to frag- 
mentation of the current layer and multiple reconnection 
sites. These instabilities have perhaps not been observed 
in previous 3D null point simulations owing to the re- 
quirement of rather large resistivity and viscosity and the 
relatively short time scales over which the boundary driv- 
ing velocities were imposed. We report here, for the first 
time, an observed instability of a Kelvin-Helmholtz (KH) 
nature occurring in a fully three-dimensional current- 
vortex sheet at the fan plane of a 3D null point. 

The KH instability has attracted significant attention 
for a number of years, both in the hy dro dy namic limit 
and in the presence of a magnetic fiel d 17 * 18 *. A uniform 
magnetic field component aligned with the shear flow is 
known to be a stabilizing influence, due to the associ- 
ated magnetic tension force. On the other hand, a mag- 
netic field component perpendicular to the shear layer (or 
'transverse field') does not affect the stability, but modi- 
fies the growth rate of a given mode^. In our simulations, 
the instability takes place at the fan plane of a 3D null, 
and there exists a strong transverse magnetic field com- 
ponent (the radial field B r associated with the potential 
field defining the null) , as well as a sheared in-plane com- 
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ponent that reverses sign at the same location as the flow 
(the azimuthal field Bq associated with the current layer 
that forms at the fan in response to the boundary driv- 
ing) . The stability of such a current- vortex sheet - where 
the magnetic and velocity shear layers coincide - has been 
previously studied in 2- and 2. 5- dimensions'^. Ein- 
audi et a/P^ showed that, in the incompressible limit, a 
transition between a tearing-like regime and a KH-like 

regime occurs when A = ) = ^' ^ ere ^ b 

and L v are the widths of the magnetic and velocity shear 
layers, Av is the velocity difference across the layer and 
ca the Alfven speed far from the layer. When A < 1, a 
tearing unstable regime is found as the magnetic shear 
strongly outweighs the velocity shear. Conversely, when 
A > 1 the velocity shear dominates the layer, and the 
linear phase of the instability is ideal, and of a KH na- 
ture. This transition has also been shown to hold in 
weakly compressible^ and viscous plasmas^. Dahlburg 
et alW* argue that even when A > 1, the presence of the 
magnetic shear fundamentally alters the nature of the 
KH-type instability by allowing magnetic reconnection 
to become important, and that therefore the instabilities 
of the current- vortex sheet should not be considered as a 
simple mix of tearing and KH modes. 

The magnetohydrodynamic (MHD) KH instability 
is important in a broad range of plasma environ- 
ments, both terrestrial and astrophysical. There are 
numerous observations of KH signatures on planetary 
magnetopauses^HZD. The KH instability is also crucial 
in the disruption of astrophysical jets, emanating for 
example from young stellar objects and active galactic 
nuclei 28 . In the solar corona there are recent direct obser- 
vations of KH instabilities in regions of strong flow shear 
associated with eruption J 29|3Q i Our work highlights that, 
in addition to these fast flows, magnetic separ at rices are 
prime locations for the formation of current- vortex sheets 
in magnetic fields of complex topology such as the so- 
lar corona. These separatrix surfaces have already been 
proposed as preferential sites of plasma heating in the 
corona^, which could be significantly enhanced by in- 
stabilities such as that studied herein. 

In this work we will investigate the self consistent for- 
mation and stability of the current-vortex sheet created 
by twisting motions around the spine foot points of a 
linear, rotationally symmetric 3D magnetic null. We re- 
strict our focus to the KH unstable regime (A > 1), leav- 
ing the tearing-type regime to a future study. It should be 
noted, however, that in the tearing unstable regime the 
growth of the instability would be expected to be slow 
as a velocity shear flow is known to damp the tearing 
mode growth rate 32 . We also note that here the geom- 
etry of the problem is much more complex than that of 
the studies discussed above. In particular, the widths of 
the shear layers {Lb and L v ) vary along the transverse 
(radial) direction, and are set in a self-consistent manner 
by a balance between the driving flow and the dissipa- 
tion in the system, rather than being fixed by the initial 



conditions. 

The investigation is structured as follows: in Section 2 
the numerical setup is discussed. In Sections 3 and 4 we 
investigate, respectively, the formation of the current- 
vortex layer, and its breakup. In Sections 5 and 6 we 
discuss our findings and present our conclusions. 



II. NUMERICAL SETUP 

The investigation was carried out numerically by solv- 
ing the compressible MHD equations in the following 
form 



V x (v x B) + r?V 2 B 
-V • (pvv) -Vp+jxB 
+ /ifv 2 v+iv(V-v) 
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with velocity v, magnetic field B, density p, thermal en- 
ergy e, gas pressure p = (7 — l)e, magnetic diffusivity 77 
and dynamic viscosity ji. The viscous heating term is 



dvi dvi 



dvj dvi 
dxi dxj 



Kv-v) ! 



using the convention of summation over repeated indices. 
In most of our simulations we use prescribed, spatially- 
uniform 77 and p. 

The equations are solved on staggered grids with a 6th 
order method applied to find the partial derivatives, a 
5th order method used for interpolation, and a 3rd or- 
der predictor-corrector method used to advance the code 
in time. In the simulations described herein, the ve- 
locity components perpendicular to the boundaries are 
set to zero, implying no mass flux through the bound- 
aries. On the x-boundaries, we restrict parallel velocities 
to be non-zero only in the two regions where a prescribed 
stressing velocity flow is imposed (see below). Near the 
y- and z- boundaries a strong damping region is in place 
(discussed below) so that in practice parallel velocities 
on these boundaries are negligible. Further information 
relating to the code can be found in Ref. [33] and on 
|http : //www . astro . ku . dk/~kg/] 

The equations are non-dimensionalised by setting the 
magnetic permeability po = 1 and the gas constant equal 
to the mean molecular weight. This results in one time 
unit representing the travel time of an Alfven wave over 
a unit distance through a plasma with unit density and 
unit magnetic field (p = 1, |B| = 1). This also means 
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that magnetic diffusivity is equal to magnetic resistiv- 
ity and takes the form of an inverse magnetic Reynolds 
number rj = = ife" 1 and the kinematic viscosity 

(y = p/p) takes the form of an inverse plasma Reynolds 
number v = = Re -1 , where L and Vq are some typ- 
ical length scale and velocity. Note that for simplicity 
we have neglected the effects of thermal conduction in 
the energy equation. This means that the temperature 
is not equilibrated along field lines. If conduction were 
included we would expect to see a modification of the 
plasma properties in the current-vortex layer that forms 
- see the discussion later. 

To investigate the properties of the current-vortex 
sheet created through twisting motions around the spine 
we start with a linear, rotationally symmetric null point 
with magnetic field B = Bo(—2x,y,z) in the center of 
a numerical box of size ±[0.25,3.5,3.5]. The grid is 
stretched to include more points near to the spine and 
fan of the null to improve the resolution of structures 
there. The plasma is assumed to be an ideal gas with 
7 = 5/3 and is initially at rest with density p = 1 and 
thermal energy e = 5/3*/ 2. Here f3* is a parameter that 
controls the plasma-/^, the ratio of plasma pressure to 
magnetic pressure: P/(B 2 /2p ) = 10p f3*/3B 2 . We set 
/3* = 0.5 and Bq = 1 in all the simulations. By set- 
ting p = 1 initially, the kinematic viscosity becomes the 
same as the dynamic viscosity (y = p/p = p) and so the 
magnetic Prandtl number may be found from Pr = p/rj. 

Because of the cylindrical symmetry of our system it is 
convenient to define a new cylindrical coordinate system 

r = r cos Oy + r sin 0z, = — r sin 9y + r cos 0z, x = x 

where r = ^ y 2 + z 2 and = tan -1 (j^j- The r and 

directions are referred to as the radial direction and 
azimuthal directions respectively. In these coordinates 
the magnetic field becomes B = rr — 2xx. Rotational 
driving on the boundaries is applied in opposite senses 
around each spine foot point of the form 

Ve (x = T0.25) = ±V (i)r (l + tanh ((1 - 36r 2 ))) (7) 

where 

V (t) = v t&nh 2 (j-j. (8) 

We choose r = 0.25 to smoothly ramp the driving up 
from zero. It has been previously noted that a near- 
discontinuous increase in driving generates fast waves, in 
addition to the main torsional Alfven wave, which focus 
on to the nuh^. We ramp the driving up slowly to avoid 
this additional complication to the evolution. 

Lastly, as we have a cylindrically symmetric system, 
we impose a cylindrical damping region beyond r = 2.8 
which removes momentum in a way which increases lin- 
early with radius (r), and limits the reflection of distur- 
bances from the y and z boundaries. Thus, if we consider 
the side boundary to be the edge of the damping region 
(r = 2.8) this boundary can be considered quasi-open. 



TABLE I. Summary of simulations 



Case 


vo 


7} 




P 


Pr — v 


/rj Resolution Stable? 


1 


0.25 5 x 10" 


-4 


1 x 10~ 4 


0.2 


160 6 


Yes 


2 


0.5 


5 x 1(T 


-4 


1 x 10~ 4 


0.2 


160 3 


Yes 


3 


1.0 


5 x 10" 


-4 


1 x 10~ 4 


0.2 


160 3 


Yes 


4 


0.5 


2 x 10- 


-4 


1 x 10~ 4 


0.5 


160 3 


Yes 


5 


0.5 


1 x 10- 


-3 


1 x 10~ 4 


0.1 


160 3 


Yes 


6 


0.5 


5 x 10" 


-4 


1 x 10~ 3 


2.0 


160 3 


Yes 


7 


0.5 


5 x 10- 


-4 


5 x 10~ 4 


1.0 


160 3 


Yes 


8 


0.5 


2 x 10- 


-4 


1 x 10~ 5 


0.05 


320 3 


No 


9 


0.5 


5 x 10- 


-4 


1 x 10~ 5 


0.02 


320 3 


No 


10 


0.5 


1 x 10- 


-3 


1 x 10~ 5 


0.01 


320 3 


No 


11 


0.5 


2 x 10- 


-4 


numerical 




320 3 


No 


12 


0.5 


5 x 10- 


-4 


numerical 




320 3 


No 


13 


0.5 


1 x 10- 


-3 


numerical 




320 3 


No 




FIG. 1. Plots viewed at t — 2.4 in the z — plane for case 3. 
The null is at the origin, the spine along y = and the fan 
at x = 0. (a) The velocity out of the plane (note the counter 
flow regions near the null), (b) the Lorentz force out of the 
plane and (c) the current density. 



III. FORMATION OF THE CURRENT-VORTEX SHEET 

A. Qualitative Behavior 

We first discuss the formation of the current-vortex 
sheet, and its structure. Several previous investiga- 
tions have observed the formation of a current layer fo- 
cused on the fan plane in response to rotational driving 
motions^^. The early stages of the simulations pro- 
ceed as follows: once the driving begins a torsional Alfven 
wave is launched from each boundary which spreads out 
as it follows the hyperbolic shape of the field toward the 
null. The current in the wave front increases as the length 
scale (perpendicular to the fan) decreases. Once both 
wave fronts get close to the fan the current diffuses into 
the fan plane itself creating a strong current layer. Of 
the studies cited above, only Galsgaard et al 34 main- 
tained the driving for long enough to see the appear- 
ance of counter rotating (i.e. against the direction of the 
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FIG. 2. Vorticity density (a) and current density (b) in the 
fan plane (x = 0) at t = 2.4 for case 3. 




FIG. 3. Sketch showing the different magnetic tension force 
in different regions close to the spine and close to the fan. 
Between the current concentrations at the fan and spine, the 
field lines pass through an ideal region where they are effec- 
tively anchored in the flow (black dots). In the two non-ideal 
regions the magnetic field lines try to straighten (red arrows) 
due to the tension force. 



driver) flow regions near the null. However, the focus 
of that investigation was on the wave dynamics so this 
was not deeply investigated and no physical reason was 
put forward for the appearance of these flows. As these 
counter flows become important for the stability of the 
current- vortex sheet, we now investigate their properties. 



B. The Counter flow 

To investigate these flows a series of simulations were 
performed at a resolution of 160 3 for various plasma and 
boundary driving parameters to observe the early dy- 
namics of the system (see table |lj cases 1 to 7 for details). 
The viscosity in these cases is chosen to be large enough 
that the current- vortex sheet remains stable, allowing us 
to focus on the formation of the counter flow regions. 
Figure [I] shows an example of the counter rotational flow 
that begins to form around t = 1.8 for case 3 (see table |l|. 
These counter rotating regions locally inhibit the veloc- 
ity shear across the fan plane and produce what we term 
a 'hole' in the vorticity layer around the null (Fig. [5|a)). 
Despite the drop in velocity shear in this region there still 
exists a strong shear in the magnetic field across the fan 



as can be seen by the strength of the current (Fig. [2^b)). 

The formation of the vortex hole arises from the in- 
terplay of forces in the vicinity of the null. Through 
symmetry considerations the magnetic tension must be 
the dominant force driving the formation of these flows. 
Comparing the Lorentz force (identical to the magnetic 
tension in the observed plane) and velocity plots in Fig. 
[I] a clear correlation between the counter flows and the 
tension force associated with the two regions of strong 
current along the spine and fan (Fig. [IJc)) is indeed ev- 
ident. The reason for the development of such a pattern 
in the tension force is not clear but is likely linked to the 
manner in which field line slippage occurs within the two 
current regions. Certainly, as the tension force is opposite 
in each current concentration, it appears that a typical 
magnetic field line can be considered as being effectively 
anchored in the ideal region between the current concen- 
trations, with this point acting as a pivot about which 
the field line wants to straighten - see Fig. [3J This sug- 
gests that the field lines around the spine have become 
partially detached from the driving boundary as a result 
of the connection change within the strong current con- 
centration near the spine. Were these foot points still 
strongly attached to the field in the volume, then the 
magnetic field would be expected to straighten around 
the foot points themselves (giving a unidirectional ten- 
sion force between the foot points and the fan plane). In 
any case, it is clear that the current concentration along 
the spine has a profound effect on the plasma dynamics 
near the fan plane. 

We find that with stronger driving (cases 1 to 3) or for 
lower values of resistivity (cases 2, 4 and 5) the current 
build-up near the spine (and therefore the magnetic ten- 
sion force in this region) increases. Thus, the strength of 
counter flows, and therefore the size of the 'hole' in the 
vorticity, also increase. In addition, when r] is reduced 
the thinning of the fan current layer combined with the 
hyperbolic shape of the magnetic field also widen the vor- 
ticity hole. As we might expect, a reduction in viscosity 
(cases 6 and 7) also increases the strength of the counter 
flows by reducing drag between the shear layers. Lastly, 
we note that the strength of the current regions near the 
driving boundaries are strongly dependent upon the cho- 
sen driving profile and that other choices of rotational 
driving will result in different degrees of counter flow. As 
we shall see later, the counter flow region plays a key role 
in determining the initial region of instability. 



IV. KH INSTABILITY OF THE CURRENT-VORTEX 
SHEET 

A. Properties of the shear layer 

Having discussed the formation of the current-vortex 
sheet at the fan surface, we proceed to explore its stability 
to the KH-type instability. Two sets of simulations were 
performed at 320 3 grid resolution, see table [l[ For one set 
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FIG. 4. Azimuthally symmetric shear layer quantities for 
case 9 at t = 5. Solid line: the fast mode Mach number 
dashed line: the projected Alfven 
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(cases 8 to 10) \i was set to zero and therefore viscosity is 
handled through numerical diffusion. This gives the least 
damping of any KH fluctuations for the chosen resolution 
and provides a benchmark against the second set (cases 
11 to 13) where /i = 1 x 1(T 5 . 

Within the fan plane current sheet the velocity shear 
is associated with a variation of vq with x. Relative to 
this flow the magnetic field has two components: a strong 
transverse guide field component B r associated with the 
initial potential null point field, and a sheared compo- 
nent Be parallel to the plasma flow (associated with the 
current in the layer formed in response to the bound- 
ary driving). Related to these two components we can 
define two radially varying but rotationally symmetric 
Mach numbers: 



M f 



Av 



Av-sjpno 

M A ,proj. = 2\£> ' 



where Av and AB are the total velocity and magnetic 
shear across the layer and c s and c A are the sound and 
Alfven speeds, all of these quantities being dependent on 
r. Mf is the fast mode Mach number related to the ve- 
locity shear and Ma, proj. is the projected Alfven Mach 
number associated with the sheared magnetic field com- 
ponent. For the KH instability, in the case of a constant 
perpendicular guide field, the instability is linearly sta- 
bilized at all wavenumbers when Mf > 2 (Ref. [T8|) . In 
addition, the transition from a tearing-type to a KH-type 

1 



2/3 

A, proj. 



instability should occur around A = (j^j M\ 

(Ref. 19). Therefore, to excite the KH instability we 
choose the driving and plasma parameters so that the 
layer is super-Alfvenic (in a projected sense) but sub- 
magentosonic. Figure [4] shows that these conditions are 
satisfied at almost all values of r for case 9 just before 
any instability arises. We note that the reduction in A 
beyond r = 2.0 is in part due to the boundary damping 
at r = 2.8, however the instability begins well away from 
this edge and so the boundaries do not affect the main 
evolution. 



B. Nature and evolution of the sheet breakup 

We now describe the growth and evolution of the KH 
instability in our current- vortex sheet. In this section we 
discuss the results of simulation 9, which are represen- 
tative of the evolution of the instability in the various 
simulations. Loosely speaking, the instability involves 
a filament ation of the current layer, with the formation 
of vortical flows in the plane of the velocity shear. Fig- 
ure [5|a-c) shows the development of the KH vortices by 
plotting the associated component of velocity perpendic- 
ular to the x = plane. In Fig. |6ja) this is quantified 
by plotting the azimuthally averaged value of pv^/2 at 
x = (the original position of the fan plane). From both 
figures it is evident that the instability initially develops 
at around r = 0.9. The vortex tubes form aligned to the 
radial magnetic field, thus reducing the damping effects 
of magnetic tension. Analyzing Fig. [6^b) it is clear that 
r = 0.9 coincides with the radial peak in vorticity. Thus, 
the size of the vortex hole in the fan plane (discussed 
above) dictates the starting point of the KH instability 
growth. We note that the development of the instability 
depletes the net vorticity within the current-vortex layer 
at later times. 

Following the guiding influence of the radial magnetic 
field the vortex tubes spread outwards and inwards from 
the initial radius of formation. As they spread outward 
from the null they encounter a longer (in the direc- 
tion) and thinner (in the x direction) shear layer and 
so branch off in order to maintain a diameter approx- 
imately equal to the thickness of the shear layer (Fig. 
|5jc)). Conversely, as they spread inwards toward the 
null they coalesce. Once out of the linear growth phase 
the KH vortices saturate as they reach the width of the 
shear layer and a new slowly varying state is reached. 

Figure [7] shows how the development of the instability 
affects the structure of the current sheet. The instability 
results in a rippling or kinking of the shear layer, in the 
same way as described in Ref. [2Ql In the strongly sheared 
stagnation point flow between each vortex a strong cur- 
rent layer forms. This fragments the current sheet into 
filaments that lie between each of the branched off vor- 
tex tubes and appear as fingers of current in Fig. [5jd-f). 
This additional localization of the current layer in the 
azimuthal direction naturally leads to the formation of 
twisted magnetic field structures (though Ampere's law) 
along each current filament (Fig. [7m) . The circular com- 
ponent of these fields is, however, small in comparison 
with the radial guide field resulting in only a small de- 
formation of the fan plane (black dashed line, Fig. [7]). 
Similar circular magnetic field structures have been ob- 
served to form between the plasma vortices in 2D simu- 
latio ns wi th comparatively better resolution of the shear 
layer^^. However, in the 2D scenario the magnetic ten- 
sion of the field in the flow vortex leads to additional cir- 
cular magnetic fields co-aligned with the vortices which 
we do not see here, perhaps as a result of the strong 
transverse field that is present. 
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FIG. 5. (a-c) The velocity out 
of the plane (v x (x = 0)). (d- 
f) The current density in the 
plane (\J\(x = 0)). The con- 
tours are scaled to the maxi- 
mum in each frame. For times 
t = 7.0 (a,d), 8.0 (b,e) and 9.0 
(c,f), for simulation 9. 




FIG. 6. Azimuthally averaged quantities plotted as a function r in the x — plane (the initial position of the fan plane) at 
different t. (a) The average perpendicular kinetic energy (< pv 2 x {x — 0) >), (b) the average vorticity density (< |w|(x = 0) >) 
and (c) the average current density (< |J|(x = 0) >), where < .. > denotes an average over the azimuthal angle = tan -1 (z/y). 
Solid line: t = 2.0, dotted: t = 3.5, dashed: t = 5.0, dot-dashed: t = 6.5, triple-dot-dashed: t = 8.0, long dash: t = 9.5. 



As a result of the formation of the flow vortices, the 
kinetic energy within the volume is more efficiently con- 
verted to local twist within the shear layer and dissipated 
through ohmic heating than prior to the onset of the in- 
stability. Figure [8] shows this for case 9 as a drop in the 
volumetric kinetic energy and viscous dissipation along 
with an increase in ohmic dissipation during the non- 
linear phase (beginning around t = 7.5). Considering 
the total dissipation within the volume it is clear that 
the onset of the instability has led to an increase in the 
localized heating of the plasma around the fan plane. 

Once the instability saturates the layer settles down to- 
ward a new slowly varying state, in which no secondary 
instabilities of the current filaments arise but some fila- 
ments are seen to begin to coalesce through the 'zipping 
up' of adjacent current branches. It was conjectured by 
Dahlburg & Einaudi 22 and subsequently confirmed by 
Onofri et alW* that a strong guide field stabilizes sec- 
ondary kinking instabilities leaving a state in which the 
coalescence instability is the dominant mode, and our re- 
sults imply that this is still the case when there is a strong 



shear flow present. Comparing the averaged shear layer 
widths and strengths from before and after the onset of 
the instability; the layer has become widened and the 
shear reduced. Thus, the KH instability acts as a key 
to the relaxation (on average) of the stress across the 
fan plane by allowing the system to transition to a state 
with a steady sharing of flux between the two topologi- 
cally distinct regions, as discussed below. 



C. Reconnection Rate 

Before going into the details of how to measure the 
reconnection rate in the layer, it is instructive to con- 
sider whether such local variations in connectivity could 
have any effect on the global connectivity of the field 
around the null. Figure [9] shows how connections near 
the null change through the onset of the KH instability 
in the layer. Initially the blue field lines (plotted from 
foot points on the driving boundaries) are connected to 
the yellow and black field lines (plotted from the foot 
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simulation 9. 
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points on the side boundaries). Once the driving begins 
the foot points of the blue field lines are advected and, 
after the main current sheet forms at the fan, these field 
lines begin to slip on to their nearest neighbors in a cir- 
cular manner. This slippage is the well studied torsional 
fan reconnect ion scenario, during which no flux is trans- 
ferred across the fan plane. The reconnection rate for 
this slippage is found from the maximum of the integral 
of E\\ along field lines passing very near to the spine and 
out along the fan. 

However, once the instability occurs in the current 
sheet numerous connection change events occur across 
the fan plane - see Fig. [9^c). That is, flux has been ex- 
changed between the two topologically distinct regions. 
How do we evaluate and interpret the reconnection rate 
in this case? Clearly, once the KH instability begins, 
both kinds of reconnection are taking place. Focusing on 
the connection change across the fan plane it is evident 
from Fig. [7] that for each KH vortex the flux transferred 
in one direction across the fan is matched by an equal 



amount transferred in the other direction, so that the 
net flux transfer is zero. This can be expressed math- 
ematically through the fact that the integral of v x B 
around a closed curve (C) in the fan plane beyond the 
edge of the current sheet (i.e. in the ideal region in which 
E = —v x B) is zero: 



v x B • dl 



E-dl 



V x E • dS 



dB 

~dt 



■ ndS = 0, 



(10) 



where S with normal vector n is the portion of the fan 
surface bounded by C, and n is perpendicular to B by 
definition. The gross rate at which flux is driven across 
the fan plane by all of the KH vortices gives the total 
associated rate of reconnection. To find this we sum the 
contribution from each half of each vortex tube in turn. 



We exploit the path independence of E (see Eq. ( 10 )) and 
consider the path depicted in Fig. [To] Due to the path 
independence, the rate at which flux is driven across the 
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FIG. 9. A selection of field lines traced from the boundaries of the box (a) initially (t = 0), (b) following the main current sheet 
formation (t = 5.0) and (c) after the onset of the instability (t = 8). The red field lines depict the position spine and fan of the 
null and the circular arrows the direction of boundary driving. A description of the connectivity change is given in the text. 





FIG. 10. Sketch of the paths considered in Eqs. (11 12). The 



grey isosurface shows the perturbed fan plane and the black 
streamlines the plasma flow. The null point is located at B. 



fan by half of the depicted vortex tube can be written as 
- / vxBd\=[ E -d\ = [ E • dl (11) 

J AC J AC J ABC 

where AB and BC are field lines lying in the fan surface, 
and so the rate of change of flux associated with half of 
one KH fluctuation, \I> kh-> is given by 



KH 



Ewdl 



AB 



Ewdl 



CB 



(12) 



Therefore, the reconnect ion rate associated with one vor- 
tex tube is double the difference between the integral of 
E\\ along the strong current lane between two adjacent 
tubes (AB) and the integral of En through the weak cur- 
rent along the tube center (CB). For the total gross rate 
of flux transfer across the fan plane this must then be 
summed over all of the vortex tubes. 

To evaluate the total reconnection rate in practice, for 
each snapshot in time, we integrate E\\ along a large num- 
ber (1800 ~ 3600) of field lines in the fan plane. Each 
integral is evaluated between r = 0.05 (since field trac- 
ing at the null is numerically problematic) and r = 2.8 



(the edge of the boundary damping region). The data 
is then binned to remove long wavelength modes result- 
ing from the Cartesian grid, after which the difference 
between peaks and troughs is summed over the entire 
angular distribution. This gives an approximation to the 
total rate of flux transfer across the fan plane 



Ri 



total 



(13) 



Consider now the rate of rotational slippage occurring 
in the volume associated with the torsional fan reconnec- 
tion that occurs even in the absence of the instability. To 
evaluate this we take advantage of the cylindrical sym- 
metry of the driving. We know that despite the small 
scale fluctuations in the main fan plane current sheet the 
driven foot points are still smoothly rotating and thus, on 
average, the flux must change connections in a symmet- 
ric manner to preserve the overall symmetry. Therefore, 
we define the rate of rotational slippage of flux as the 
maximum of the azimuthal average of the integral of E\\, 
i.e. 



((/ 



Endl 



(14) 



where < .. > denotes an average over the azimuthal an- 
gle. In practice this occurs along field lines lying asymp- 
totically close to the spine and fan. Note that this defines 
an overall reconnection rate, including rotational slippage 
both within the fan current sheet, and the large-scale cur- 
rent concentrations around the spine. Using these defini- 
tions we can quantify how flux is reconnected as various 
parameters of the system are varied. 

In general, in each simulation, \I/ S grows steadily un- 
der the action of the continued boundary driving until, at 
some point, the KH instability fragments the the current- 
vortex layer, decreasing \I/ S as Rt grows rapidly. Fig- 
ure ll'b-c) shows this for cases 8-10 where it is clear 
that Rt rapidly grows to dominate over \I/ S as a result 
of the recursive nature of the KH induced connection 
change and the large number of reconnection sites in the 
fan plane. 
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FIG. 12. Scalings with 77 of the (a) average shear layer width 
(L v , between r = 0.6 and 2.8), (b) radius of initial instability 
growth (r inst .), (c,d) growth rates of (K x , R T ) T Rt ), 
normalized by L/ca), (e,f) growth rates (K x , Rt) (normalized 
by (L/ Av)min), (g,h) peak values of (K x , Rt)- Asterisks: 
with numerical viscosity, diamonds: ji = 1 X 10 -5 . 



D. Quantitative properties of the system 

Here we explore briefly how the behavior described 
above changes as we vary the resistivity and viscosity of 
the plasma. As the resistivity of astrophysical plasmas 
such as the solar corona is incredibly small (and there- 
fore impossible to simulate numerically) it is important 
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FIG. 13. Current density in the x = plane for (a) case 13 
(at t = 6.5) and (b) case 11 (at t = 7.5). 



that we understand how any simulated reconnecting sys- 
tem behaves as the resistivity is reduced. Typically, this 
involves finding scaling laws for parameters such as the 
growth rate of the instability or the peak reconnection 
rate. A full scaling analysis is impractical here given the 
size of the data sets and computational power required, 
but we have made a preliminary investigation. 



1. Quantitative properties of the current- vortex sheet 

We first consider how the properties of the current- 
vortex sheet changes as we vary rj. As r] is reduced the 
current layer, and thus the velocity shear layer, becomes 
thinner. This can be seen in Fig. [l2^a) which shows the 
mean velocity layer thickness (L v ) between 0.6 < r < 2.8 
in cases 8 to 13 at t = 4. As the shear layer thins, 
the vorticity hole becomes wider (due to the hyperbolic 
shape of the field) and the KH instability sets in initially 
at a larger radius (Fig. [l2|fo)). This suggests that at 
realistic coronal parameters the instability of the sheet 
will occur at large distances from the null, but through 
the spreading of its influence described in the previous 
section will still dominate a significant portion of the fan 
plane current sheet. In competition with the above effect, 
the increase in viscosity between cases 8 to 10 and 11 to 
13 clearly widens the shear layer so that in more viscous 
fluids the vortices form closer to the null. 

An important consequence of the formation of a thin- 
ner shear layer is that the number of vortices that form 
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greatly increases (Fig. 13). This suggests that as 77 is re- 
duced and the layer thins the associated gross flux trans- 
fer across the fan plane may increase. We will explore 
this further below. Lastly, at large values of 77 the mag- 
netic field may slip through the plasma more readily than 
at lower values. Thus, as 77 is lowered, so the relative ve- 
locity to magnetic shear (M^, pr oj.) decreases. This is 
shown in Fig. [14] which also shows the increasing current 
and vorticity density as a result of the thinning of the 
layer. 



2. Quantitative properties of the instability 

We now consider the dependence of the instability on 
77 and v. Note that for a KH-type instability we would 
expect no dependence on 77 in the linear phase. However, 
there is an indirect influence due to the varying prop- 
erties of the shear layer that forms, as discussed above. 
Moreover, the properties of the instability as it saturates 
can be expected to depend on the dissipation parameters, 
as discussed below. 

The evolution and growth of the instability of the 
current-vortex layer can be quantified by integrating the 
kinetic energy associated with the component of the flow 
perpendicular to the x = plane, i.e. 



K x 



pv^dydz, 



(15) 



evaluated in the x = plane. This can be compared with 
the gross rate of flux transfer across the fan surface (Rt- 
Eq. 13). As shown in Fig. 11, both quantities exhibit an 



exponential growth phase during the early stages of the 
instability. The increase in Rt lags behind the growth of 
K x as expected for a KH-type instability since the linear 
phase of the instability is ideal, and initially the kinetic 
energy of each vortex is expended in ideally deforming 
the fan surface (dashed black line, Fig. u\. By assuming 
that the growth phase approximates ~ e 1 ^ and normal- 
izing against the Alfven travel time across the width of 
the simulation box (using the Alfven speed at the spine 
foot points and the distance between the ^-boundaries: 
tA ~ L/ca(x = 0.25) = 0.5/0.5 = 1) we can compare 
the growth rates of K x and Rt (denoted V and Tr t , re- 
spectively) relative to a typical time scale for the whole 
system. It should be noted, however, that as these growth 



rates are global quantities they inherently include many 
different effects that can affect the growth rate of the 
KH instability (such as different dominant wave numbers, 
shear layer widths, Alfven Mach numbers and magnetic 
and plasma Reynolds numbers of each shear layer). 

Fig. [T2fc) shows the growth rate of K x for the different 
cases that were studied. In general, the growth rate of 
the instability is seen to increase as 77 is reduced. The 
exception to this is case 8 (with \i = 1 x 10 -5 and 77 = 
2 x 10 -4 ), where the wavelength of the dominant wave- 
mode has likely become comparable to the viscous length 
scale. The growth rates of Rt (Fig. 12 ^d)) follow a 
similar trend, but at the lowest values of 77 are reduced 
as a likely result of increased numerical diffusion leading 
to an underestimation of the value of 77. 

One of the factors affecting the growth of the instabil- 
ity that we can normalize against is the changing thick- 
ness of the shear layer. Specifically, we can consider the 
effect of normalizing against the minimum travel time 
across the vorticity layer ((L v / Av) m i n ), taken before the 
onset of the instability (at t = 4). Fig. [l2^e,f) shows the 
same growth rates as above but with the effect of the 
thinning shear layer removed. In this case the growth 
rates of both quantities now decrease with decreasing 77. 
This shows that the dominant factor in the increase in 
the global growth of the KH instability as 77 is reduced 
is the thinning of the shear layer. What is clear from 
our results is that the growth rates of both K x and Rt 
are linked and that the relative strength of resistivity to 
viscosity (i.e. the magnetic Prandtl number) plays an 
important role in setting up the initial layer, and thus 
how fast the KH instability of the layer grows. 

Lastly, we note that small-scale near-turbulent recon- 
nection events around the edge of the flow vortices have 
been seen in 2D simulations (e.g. Ref. 21) to help drive 
the plasma circulation and increase the growth rate. This 
effect is beyond what is currently available to study at 
the resolution of these simulations but could also play 
an important role in the dynamics of the current- vortex 
sheet. 

Once the vortices grow to be comparable to the width 
of the shear layer they saturate and reduce in amplitude. 
At this point both K x and Rt peak, before reaching a 
new steady level. Fig. [l2|g,h) shows the peak values at- 
tained. There appears to be a reduction in the peak value 
of K x with decreasing 77, although this does not appear 
systematic and may be partly due to an under-resolution 
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of the vortex structures in v x at the smallest values of rj. 
It could also be that the maximum energy content of the 
perturbation flow associated with the vortices is reduced 
as the non-linear phase is entered earlier with a reduction 
in the width of the shear layer. By contrast, the peak re- 
connection rate Rt increases approximately linearly as rj 
is reduced. This surprising result arises because as r] is 
reduced there are many extra vortex tubes produced as 
the shear layer thins. 



V. DISCUSSION 

In this work we investigated first the self consistent 
formation, and then the stability to Kelvin-Helmholtz- 
type vortices, of the current- vortex sheet which forms 
in response to twisting motions around the spines of a 
symmetric 3D magnetic null point. The formation of the 
current-vortex layer was found to be complex, with local 
non-ideal effects brought about by strong current regions 
that form along the spine lines. This strongly affects 
the initial region of KH growth when the layer is KH 
unstable. 

The instability leads to the breakup of the planar cur- 
rent layer at the fan, into vortex structures in the plane 
of the flow and magnetic shears, consistent with earlier 
studies of 2- and 2.5-dimensional current- vortex sheets. 
Due to a combination of the field geometry and resis- 
tive and viscous diffusion, the unstable region with high 
magnetic and flow shears is confined at intermediate radii 
from the null point. The widths and magnitudes of the 
current and vorticity shear layers are consistent with a 
KH-type instability according to the previous theor}D^. 
There are a number of aspects of our simulations that 
point to the observed instability being of a KH-type. 
First, the growth of the reconnected flux (Rt) lags be- 
hind the growth of the kinetic energy of the perturbation, 
suggesting that in the early stages the instability is pre- 
dominantly an ideal one. Furthermore, the instability is 
associated with a rippling or kinking of the shear layer, 
consistent with the KH-dominated regime in 2D as dis- 
cussed in Ref. [20l One new effect that we observed was a 
'branching' of the vortex tubes / current filaments in the 
direction transverse to the shear. This is an effect that 
arises due to the fully 3D nature of the system, in that 
the layer gets longer (in the azimuthal direction of the 
shear flow) and thinner as we move away from the null 
in this transverse direction. 

Once formed, it was found that the instability quickly 
grows to dominate the majority of the current-vortex 
sheet. The strong current lanes between the 'branched' 
vortex tubes more efficiently heat the plasma in the layer 
and enable a rapid recursive sharing of magnetic con- 
nections between the two previously separate topological 
regions. This associated connectivity change dominates 
over the torsional slippage associated with the global 
driving flow. This shows how, subject to a smooth de- 
formation, instabilities in the fan surface current layer 



of 3D null points can lead to a sudden burst of mix- 
ing of plasma between the previously distinct regions on 
either side of the separatrix surface. Such recursive re- 
connection is also a feature of reconnection between nulls 
connected by multiple separators^. Plasma heating at 
fan separatrix surfaces has previously been proposed to 
be important in the solar corona 3 -^. Our results suggest 
that this heating could be significantly enhanced by in- 
stabilities of the current vortex-sheet that forms at the 
separatrix. It should be noted that the effect of ther- 
mal conduction in these simulations has been neglected. 
This would act to transport thermal energy away from 
the current layer along field lines, and thus reduce the 
temperature, altering the pressure/density distribution 
within the current- vortex layer. While we do not expect 
that our qualitative results would change, we might ex- 
pect the region of the fan plane that is unstable, and the 
growth rate, to be slightly modified. The effect of addi- 
tional terms in the energy equation should be considered 
in future studies. 



As the instability of the current-vortex sheet is strongly 
dependent on the driving and plasma parameters we have 
only focused on a small parameter range to best optimize 
the available resolution. We focused our attention on the 
importance of viscosity and resistivity in the dynamics of 
the layer. We were hampered by the difficulty in resolv- 
ing the fine scale structure following the formation of the 
layer. However, what became clear is that both the resis- 
tivity and the viscosity of the plasma have a strong effect 
on the initial position, growth rate and saturation level of 
the KH instability in the fan plane current- vortex sheet. 
This is perhaps counter-intuitive since the KH instability 
is an ideal one, but the key point is that the properties of 
the shear layer itself are strongly dependent on 77 and z/, 
which means that the growth rate of the instability also 
is, albeit indirectly. In addition, changing other parame- 
ters would also affect the layer. For instance in 2D it is 
known that increasing the compressibility of the plasma 
or the driving speed so that Mf > 1.3 produ ces sh ocks 
where the plasma is accelerated by the vortices^^J 



We note that for continued driving beyond the end of 
these simulations the kink instability may set in around 
the spine lines which would destroy the formation of the 
current- vortex shear layer, leading to patchy reconnec- 
tion across the fan plane 3 . We would expect that the 
kink instability sets in sooner for systems which are more 
strongly driven or less resistive, as more twist may build 
up around the spine. This suggests that there may only 
be a window in time for which the KH instability forms 
in the fan plane current- vortex sheet before the kink in- 
stability destroys the shear layer. It remains to be seen 
if at more realistic coronal parameters this window is too 
short for the KH instability, under these driving condi- 
tions, to be a realistic means for sudden energy release. 
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VI. CONCLUSIONS 

In future there are a number of extensions to our study 
that could be considered. Relaxing the symmetry of ei- 
ther the initial magnetic null point field or the driving 
flow is likely to profoundly affect the locations in which 
the current- vortex sheet is unstable. In addition, we have 
so far not probed the range of parameters appropriate to 
a tearing- type instability. Furthermore, we have observed 
a similar type of instability when the null is subjected 
to shearing rather than the highly-symmetric rotational 
perturbations considered here. The resulting spine-fan 
reconnection mode is arguably more generic, and thus it 
may be a good candidate for rapid energy release through 
either the tearing or KH instabilities. This will form the 
basis of a future study. The development of a more tur- 
bulent state should also be studied, though this requires 
much higher resolution of the current-vortex sheet, and a 
different numerical approach may be necessary. The in- 
stability of current layers around 3D null points will also 
have a profound effect on the acceleration of particles in 
these locations, with the well-collimated beams observed 
in existing studied likely to be modified. 

We conclude that the fragmentation via the KH in- 
stability of the torsional fan current-vortex sheet could 
provide a rapid mechanism for energy release and con- 
nectivity change between the two topologically distinct 
regions separated by the fan separatrix surface. How- 
ever, the setup considered here is highly symmetric and 
other instabilities may be triggered first if wider param- 
eter spaces are considered. Further work is needed to 
determine if this mechanism is a realistic candidate for 
sudden energy release in parameter spaces typical of as- 
trophysical plasmas. We hope that the techniques and 
understanding developed in this work could be used to 
probe this scenario more deeply in the future. 
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